Motif / TFBS Analysis

Workflow to ID motifs and match to genes

(Full implementation in Snakefile)

1. Download TF motifs from cisbp
2. Only keep TF which are abs(logFc) > 1 between GFP/RPE and iPSC (~1200)
3. Check for motifs in narrow ATAC-seq peaks with fimo
4. Identify closest 2 genes (under 500k bp) to each motif
5. Bootstrap steps 2 and 3 250 times each to get background rate

# Load Libraries without printing any warnings or messages
library(tidyverse)
library(ggridges)
library(cowplot)
library(DESeq2)
library(ggrepel)
load('/Volumes/data/projects/nei/hufnagel/iPSC_RPE_ATAC_Seq/Rdata/dds.Rdata')
load('/Volumes/data/projects/nei/hufnagel/iPSC_RPE_ATAC_Seq/Rdata/sample_bootstrap_counts.Rdata')
load('/Volumes/data/projects/nei/hufnagel/iPSC_RPE_ATAC_Seq/Rdata/sample_closestTSS.Rdata')

DESeq2-based analysis

I realized that motifs by samples is similar to genes by samples

GFP vs iPSC

dds$Type <- relevel(dds$Type, ref='RFP')
res_GFP_iPSC <- results(dds, contrast=c("Type","GFP","iPSC")) 
res_GFP_RFP <- results(dds, contrast=c("Type","GFP","RFP")) 
#plotMA(res_GFP_iPSC)
#plot(hist(res_GFP_iPSC$pvalue, breaks=1000))

# Load in RNA-seq data
gfp_rfp <- read_csv('~/git/ipsc_rpe_RNA-seq/data/GFP_vs_RFP.results.csv')

Parsed with column specification:
cols(
  Gene = col_character(),
  baseMean = col_double(),
  log2FoldChange = col_double(),
  lfcSE = col_double(),
  stat = col_double(),
  pvalue = col_double(),
  padj = col_double()
)

rpe_ipsc <- read_csv('~/git/ipsc_rpe_RNA-seq/data/RPE_vs_iPSC.results.csv')

Parsed with column specification:
cols(
  Gene = col_character(),
  baseMean = col_double(),
  log2FoldChange = col_double(),
  lfcSE = col_double(),
  stat = col_double(),
  pvalue = col_double(),
  padj = col_double()
)

#gfp_rfp %>% head()

Table of results

Motifs high in GFP

# calculate z scores for sample_bootstrap_counts
Z_scoring <- sample_bootstrap_counts %>% 
  # collapse by Type, motif, and whether bootstrap or real
  group_by(sample, motif_alt_id) %>% 
  mutate(`Z score` = scale(Count)) %>% 
  filter(bootstrap == 'real') %>% 
  select(-bootstrap, -TF) %>% 
  left_join(tf_motif)

Joining, by = "motif_alt_id"

#Z_scoring %>% head()

The Z score is the number of standard deviation of motifs found in the sample peaks above the random set of peaks

results <- res_GFP_iPSC %>% 
  data.frame() %>% 
  rownames_to_column('motif_alt_id') %>% 
  left_join(tf_motif) %>% 
  arrange(pvalue) %>% 
  data.frame() %>% 
  left_join(Z_scoring %>%  mutate(Type = case_when(grepl('GFP', sample) ~ 'GFP',
                                                       grepl('RFP', sample) ~ 'RFP',
                                                       TRUE ~ 'iPSC')) %>% 
              filter(Type=='GFP') %>% group_by(motif_alt_id) %>% summarise(`GFP Z score` = mean(`Z score`))) %>% 
  left_join(Z_scoring %>%  mutate(Type = case_when(grepl('GFP', sample) ~ 'GFP',
                                                        grepl('RFP', sample) ~ 'RFP',
                                                        TRUE ~ 'iPSC')) %>% 
              filter(Type=='iPSC') %>% group_by(motif_alt_id) %>% summarise(`iPSC Z score` = mean(`Z score`))) %>% 
  left_join(gfp_rfp %>% mutate(TF=Gene, log2FC_GFP_RFP_RNASeq = log2FoldChange) %>% select(TF, log2FC_GFP_RFP_RNASeq)) 

Joining, by = "motif_alt_id"
Joining, by = "motif_alt_id"
Joining, by = "motif_alt_id"
Joining, by = "TF"

results %>% filter(`iPSC Z score` < `GFP Z score`) %>% arrange(pvalue) %>% DT::datatable()

Motifs high in iPSC

results %>% filter(`iPSC Z score` > `GFP Z score`) %>% arrange(pvalue) %>% DT::datatable()

Volcano

volcano_maker <- function(df, title){
  df$Class <- 'Not significant'
  df$Class[df$padj < 1e-100 & df$log2FoldChange > 1.5] <- "FDR < 1e-100 &\nlog2FC > 1.5"
  df$Class[df$padj < 1e-100 & df$log2FoldChange < -1] <- "FDR < 1e-100 &\nlog2FC < -1"
  df$Class <- factor(df$Class, levels=c('Not significant', "FDR < 1e-100 &\nlog2FC > 1.5", "FDR < 1e-100 &\nlog2FC < -1"))
  plot <- ggplot(data=df,aes(x=log2FoldChange,y=-log10(pvalue))) + 
    #geom_point(aes(colour=Class, size=`Z score`), alpha=0.5) +
    geom_point(aes(colour=Class), alpha = 0.5) + 
    scale_colour_manual(values=c("gray","darkred", "darkblue")) + 
    geom_text_repel(data=df %>% filter((padj < 1e-100  & log2FoldChange > 1.5) | (padj < 1e-100 & log2FoldChange < -1)), 
                    aes(label=TF)) +
    # geom_vline(aes(xintercept=-0.5),linetype="dotted") +
    # geom_vline(aes(xintercept=0.5),linetype="dotted") +
    scale_x_continuous(breaks=c(-2,-1.5,-1,-0.5,0,0.5,1,1.5,2)) +
    ggtitle(title) + theme_minimal()
  return(plot)
}
volcano_maker(results, 'GFP+ RPE vs iPSC TFBS motif counts')

Plots of the top motifs more common in GFP than iPSC

# motif ggridges plotter
# removes GFP 
plotterGi <- function(motif, scale=TRUE) {
  if (scale == TRUE){
    sample_bootstrap_counts <- sample_bootstrap_counts %>% filter(motif_alt_id == !!motif) %>% 
      filter(!grepl('RFP',sample)) %>% 
      group_by(sample, motif_alt_id) %>% 
      mutate(`Z score` = scale(Count)) %>% 
      ungroup() 
    sample_bootstrap_counts %>% 
      filter(!grepl('RFP',sample)) %>% 
      #mutate(sample = factor(sample, levels=c('iPSC_IIi_9','iPSC_IIJ_10','iPSC_IIK_11','iPSC_IIL_12','RFP_IIE_1','RFP_IIF_2','RFP_IIG_3','RFP_IIH_4','GFP_IIE_5','GFP_IIF_6','GFP_IIG_7','GFP_IIH_8'))) %>% 
      mutate(sample = factor(sample, levels=c('iPSC_IIi_9','iPSC_IIJ_10','iPSC_IIK_11','iPSC_IIL_12','GFP_IIE_5','GFP_IIF_6','GFP_IIG_7','GFP_IIH_8'))) %>% 
      ggplot(aes(x=`Z score`, y=sample)) + 
      geom_density_ridges() +
      geom_point(data = sample_bootstrap_counts %>% 
                   filter(bootstrap == 'real', motif_alt_id == !!motif) %>% 
                   ungroup(), aes(x=`Z score`,y=sample), colour='blue', size=2, alpha=0.5) +
      scale_y_discrete(expand = c(0.01, 0)) +
      theme_ridges() + 
      ggtitle(paste0(motif, ' (', 
                     tf_motif %>% filter(motif_alt_id == !!motif) %>% pull(TF),
                     ')'))
  }
  else {
    sample_bootstrap_counts %>% filter(motif_alt_id == !!motif) %>% 
      filter(!grepl('RFP',sample)) %>% 
      #mutate(sample = factor(sample, levels=c('iPSC_IIi_9','iPSC_IIJ_10','iPSC_IIK_11','iPSC_IIL_12','RFP_IIE_1','RFP_IIF_2','RFP_IIG_3','RFP_IIH_4','GFP_IIE_5','GFP_IIF_6','GFP_IIG_7','GFP_IIH_8'))) %>% 
      mutate(sample = factor(sample, levels=c('iPSC_IIi_9','iPSC_IIJ_10','iPSC_IIK_11','iPSC_IIL_12','GFP_IIE_5','GFP_IIF_6','GFP_IIG_7','GFP_IIH_8'))) %>% 
      ggplot(aes(x=Count, y=sample)) + 
      geom_density_ridges() +
      geom_point(data = sample_motifs %>% 
                   filter(motif_alt_id == !!motif) %>% 
                   group_by(sample, motif_alt_id) %>%
                   summarise(Count=n()) %>% ungroup(), aes(x=Count,y=sample), colour='blue', size=2, alpha=0.5) +
      scale_y_discrete(expand = c(0.01, 0)) +
      theme_ridges() + 
      ggtitle(paste0(motif, ' (', 
                     tf_motif %>% filter(motif_alt_id == !!motif) %>% pull(TF),
                     ')'))
  }
}

plots <- list()
for (i in results %>% filter((`GFP Z score`) > 5) %>% arrange(pvalue) %>% head(n=12) %>% pull(motif_alt_id)){
  plots[[i]] <- plotterGi(i, scale=T)
}
plot_grid(plotlist = plots, ncol=3)

Picking joint bandwidth of 0.176
Picking joint bandwidth of 0.168
Picking joint bandwidth of 0.166
Picking joint bandwidth of 0.159
Picking joint bandwidth of 0.159
Picking joint bandwidth of 0.168
Picking joint bandwidth of 0.162
Picking joint bandwidth of 0.167
Picking joint bandwidth of 0.163
Picking joint bandwidth of 0.164
Picking joint bandwidth of 0.166
Picking joint bandwidth of 0.166

Plots of the top motifs more common in iPSC than GFP

plots <- list()
for (i in results %>% filter((`iPSC Z score`) > 5) %>% arrange(pvalue) %>% head(n=6) %>% pull(motif_alt_id)){
  plots[[i]] <- plotterGi(i, scale=T)
}
plot_grid(plotlist = plots, ncol=3)

Picking joint bandwidth of 0.217
Picking joint bandwidth of 0.223
Picking joint bandwidth of 0.22
Picking joint bandwidth of 0.241
Picking joint bandwidth of 0.241
Picking joint bandwidth of 0.243

Closest TSS for OTX2 (M5700_1.02)

Two scoring systems:

1. Mean Count: number of motifs linked to gene (closest two genes within 500kb)
2. Sum -log10(p value): sum of the -log10(p value) of the specificity of the motif matching. This score weighs better matching motifs more than counts.

Filtering on log2(baseMean) > 5 (~ top2/3 of genes)

plotter('M0952_1.02')

Closest TSS for JUND (M4464_1.02)

sample_closestTSS %>% 
  #filter(fimo_pvalue < 1e-5) %>% 
  filter(motif %in% c('M4464_1.02')) %>% 
 mutate(Type = case_when(grepl('GFP', sample) ~ 'GFP',
                           grepl('RFP', sample) ~ 'RFP',
                           TRUE ~ 'iPSC')) %>% 
  group_by(Gene, sample, Type, motif) %>%
  summarise(Count=n(), motifs = list(motif_loc), scaleP = sum(-log10(fimo_pvalue))) %>%
  group_by(Gene, Type, motif) %>% 
  summarise(`Mean Count` = mean(Count), `Sum -log10(p value)` = sum(scaleP), motifs = list(unique(motifs))) %>% 
  arrange(-`Mean Count`) %>% 
  left_join(., gfp_rfp %>% mutate(log2FC_GFP_RFP_RNASeq = log2FoldChange, baseMean_GFP_RFP = baseMean) %>% select(Gene, log2FC_GFP_RFP_RNASeq, baseMean_GFP_RFP)) %>% 
  left_join(., rpe_ipsc %>% mutate(log2FC_RPE_iPSC_RNASeq = log2FoldChange, baseMean_GFP_iPSC = baseMean) %>% select(Gene, log2FC_RPE_iPSC_RNASeq, baseMean_GFP_iPSC)) %>% 
  filter(log2(baseMean_GFP_iPSC) > 5) %>% 
  select(Gene, Type, motif, `Mean Count`, `Sum -log10(p value)`, log2FC_GFP_RFP_RNASeq, log2FC_RPE_iPSC_RNASeq, motifs) %>% 
  head(2000) %>% 
  DT::datatable()

Joining, by = "Gene"
Joining, by = "Gene"
It seems your data is too big for client-side DataTables. You may consider server-side processing: https://rstudio.github.io/DT/server.htmlIt seems your data is too big for client-side DataTables. You may consider server-side processing: https://rstudio.github.io/DT/server.html

Closest TSS for SNAI1 (M6468_1.02)

sample_closestTSS %>% 
  #filter(fimo_pvalue < 1e-5) %>% 
  filter(motif %in% c('M6468_1.02')) %>% 
 mutate(Type = case_when(grepl('GFP', sample) ~ 'GFP',
                           grepl('RFP', sample) ~ 'RFP',
                           TRUE ~ 'iPSC')) %>% 
  group_by(Gene, sample, Type, motif) %>%
  summarise(Count=n(), motifs = list(motif_loc), scaleP = sum(-log10(fimo_pvalue))) %>%
  group_by(Gene, Type, motif) %>% 
  summarise(`Mean Count` = mean(Count), `Sum -log10(p value)` = sum(scaleP), motifs = list(unique(motifs))) %>% 
  arrange(-`Mean Count`) %>% 
  left_join(., gfp_rfp %>% mutate(log2FC_GFP_RFP_RNASeq = log2FoldChange, baseMean_GFP_RFP = baseMean) %>% select(Gene, log2FC_GFP_RFP_RNASeq, baseMean_GFP_RFP)) %>% 
  left_join(., rpe_ipsc %>% mutate(log2FC_RPE_iPSC_RNASeq = log2FoldChange, baseMean_GFP_iPSC = baseMean) %>% select(Gene, log2FC_RPE_iPSC_RNASeq, baseMean_GFP_iPSC)) %>% 
  filter(log2(baseMean_GFP_iPSC) > 5) %>% 
  select(Gene, Type, motif, `Mean Count`, `Sum -log10(p value)`, log2FC_GFP_RFP_RNASeq, log2FC_RPE_iPSC_RNASeq, motifs) %>% 
  head(2000) %>% 
  DT::datatable()

Joining, by = "Gene"
Joining, by = "Gene"
It seems your data is too big for client-side DataTables. You may consider server-side processing: https://rstudio.github.io/DT/server.htmlIt seems your data is too big for client-side DataTables. You may consider server-side processing: https://rstudio.github.io/DT/server.html

Gene-centered approach

What motifs are associated with a particular gene?

ABCA4 for this example

sample_closestTSS %>%
  mutate(Type = case_when(grepl('GFP', sample) ~ 'GFP',
                          grepl('RFP', sample) ~ 'RFP',
                          TRUE ~ 'iPSC')) %>% 
  filter(Gene == 'ABCA4') %>% 
  group_by(Gene, sample, Type, motif) %>%
  summarise(Count=n(), motifs = list(motif_loc), scaleP = sum(-log10(fimo_pvalue))) %>% 
  group_by(Gene, Type, motif) %>% 
  summarise(`Mean Count` = mean(Count), `Sum -log10(p value)` = sum(scaleP), motifs = list(unique(motifs))) %>% 
  arrange(-`Mean Count`) %>% 
  left_join(gfp_rfp %>% mutate(log2FC_GFP_RFP_RNASeq = log2FoldChange) %>% select(Gene, log2FC_GFP_RFP_RNASeq)) %>% 
  left_join(rpe_ipsc %>% mutate(log2FC_RPE_iPSC_RNASeq = log2FoldChange) %>% select(Gene, log2FC_RPE_iPSC_RNASeq)) %>% 
  left_join(., tf_motif %>% mutate(motif = motif_alt_id)) %>% 
  select(Gene, Type, motif, TF, `Mean Count`, `Sum -log10(p value)`, log2FC_GFP_RFP_RNASeq, log2FC_RPE_iPSC_RNASeq, motifs) %>%   
  DT::datatable()

Joining, by = "Gene"
Joining, by = "Gene"
Joining, by = "motif"

# library(tidygraph)
# library(ggraph)
# 
# node_data <- sample_closestTSS %>% 
#  # filter(fimo_pvalue < 1e-6) %>% 
#   #filter(motif == 'M5700_1.02') %>% 
#   mutate(Type = case_when(grepl('GFP', sample) ~ 'GFP',
#                            grepl('RFP', sample) ~ 'RFP',
#                            TRUE ~ 'iPSC')) %>% 
#   group_by(Gene, sample, Type) %>%
#   summarise(Count=n(), motifs = list(motif_loc)) %>%
#   group_by(Gene, Type) %>% 
#   summarise(`Mean Count` = mean(Count)) %>% 
#   arrange(-`Mean Count`) %>% 
#   filter(`Mean Count` > 3.9) %>% 
#   rowid_to_column("id")
# 
# 
# edge_data <- node_data %>% mutate(to = id, from = 33, weight= `Mean Count`) %>% ungroup() %>% 
#   select(from, to, weight, Type)
# 
# node_data <- bind_rows(node_data, tibble(id = 33, Gene = 'OTX2', Type = 'TF', `Mean Count` = 0))
# 
# 
# routes_tidy <- tbl_graph(nodes = node_data %>% mutate(node = as.character(id)), edges = edge_data %>% mutate(from=as.character(from), to=as.character (to)), directed = TRUE)
# 
# routes_igraph <- graph_from_data_frame(d = edge_data, vertices = node_data, directed = TRUE)
# 
# ggraph(routes_tidy) + 
#   geom_edge_link(aes(width = weights, color=as.factor(Type))) + 
#   scale_edge_width(range = c(0.2, 2)) + geom_node_point() + theme_graph() + geom_node_text(aes(label = Gene), repel = TRUE) 

# What genes have more PAX6 'associated' motifs compared from GFP to RFP
# enriched_genes <- both %>% 
#   # only keep one gene per motif
#   group_by(motif_loc, sample, Gene) %>% 
#   top_n(1, distance) %>% 
#   ungroup() %>% 
#   # keep up to two genes per motif
#   group_by(motif_loc, sample) %>% 
#   top_n(2, distance) %>% 
#   ungroup() %>% 
#   # arrange by genes most linked to motif  
#   group_by(Gene, sample) %>% 
#   summarise(Count=n(), paste(motif_loc, collapse=', ')) %>% 
#   ungroup() %>% 
#   # collapse to GFP/RFP/iPSC
#   mutate(Type = case_when(grepl('GFP', sample) ~ 'GFP',
#                           grepl('RFP', sample) ~ 'RFP',
#                           TRUE ~ 'iPSC')) %>% 
#   group_by(Gene, Type) %>% 
#   summarise(Total=sum(Count)) %>% 
#   arrange(-Total) %>% 
#   spread(Gene, Total) %>% t() 
# 
# colnames(enriched_genes) <- enriched_genes[1,]
# enriched_genes <- enriched_genes[-1,] %>% data.frame() %>% rownames_to_column('Gene') %>%  mutate(GFP = as.numeric(GFP), iPSC = as.numeric(iPSC), RFP = as.numeric((RFP)))
# 
# enriched_genes[is.na(enriched_genes)] <- 0
# 
# enriched_genes %>% mutate(`deltaGFP <-> RFP` = GFP - RFP) %>% arrange(-`deltaGFP <-> RFP`, GFP) %>% head(1000) %>% DT::datatable(rownames = F)

## does PAX6 regulate PAX6?
#Yes, yes it does.
#GFP specific!
  
# both %>% 
#   # only keep one gene per motif
#   group_by(motif_loc, sample, Gene) %>% 
#   top_n(1, distance) %>% 
#   ungroup() %>% 
#   # keep up to two genes per motif
#   group_by(motif_loc, sample) %>% 
#   top_n(2, distance) %>% 
#   ungroup() %>% 
#   # arrange by genes most linked to motif  
#   group_by(Gene, sample) %>% 
#   summarise(Count=n(), `Motif Locations` = paste(motif_loc, collapse=', ')) %>% 
#   arrange(-Count) %>% 
#   filter(Gene=='PAX6')